Student Name(s): Pratyush Kumar
Student Number(s): 5359252
In this assignment we do a time series analysis of NDVI using the HANTS (Harmonic Analysis of Time Series) algorithm and visualize the results in Python and QGIS.
For a description of the assignment consult the CIE4604-A3-NDVI_HANTS.pdf document.
In your report (or here) include at least the following results:
Part I.a
Part I.b
Furthermore, include in the zip file accompanying the report,
For extraction of the raster images by the extent defined by the shapefile, i did the following:
clip raster by mask layerraster calculator was used to produce the NDVI. The colormap used for this purpose, is a diverging colormap with shades of brown for values below 0 and shades of green for values above 0. This was chosen also in a color blind freindly version of brown and green so that the maps are visible to everyone as intended. Moreover, the color scheme was chosen such that the greener pixel is indicative of higher NDVI value.
The GDAL command to make raster masks were done as follows

The image below is of the layout formed using QGIS:

The image below is of the attribute table exported as an xlsx file:

HANTS is a software using the Fourier Analysis in order to detect outliers (clouds) and reconstruct image time series.
The analysis is performed over Sector BXII, an agricultural area located close to Sevilla (Spain).
Input data files for this script are:
NDVI_BXII_time_series_2017_2018_S2A_SP.tif (Geotiff with NDVI time series)NDVI_BXII_time_series_images_list.csv (Csv-file with the date, day number and day of year for the images)Optional (non-required) input file(s) are:
Training_2017_UTM30N_WGS84_SP.shp (Training dataset)The actual HANTS algorithm is implemented in the Python function hants.py, see also help(hants.py).
# Importing all the libraries required for this assignment
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colorbar as colorbar
import matplotlib.colors as colors
import matplotlib.animation as animation
import matplotlib as mpl
import pandas as pd
import rasterio
import shapefile
from scipy.io import loadmat
from datetime import date, datetime
from hants import hants
The input files are:
NDVI_BXII_time_series_2017_2018_S2A_SP.tif - geotiff with NDVI NDVI_BXII_time_series_images_list.csv (Csv-file with the date, day number and day of year for the images)The code for reading these dataset was already developed in Exercise 4. The following variables are created
You should definitely take a look at each variable to see what it is about.
# set the filenames and filepath
path = './input_files/Hants/'
NDVIfile = 'NDVI_BXII_time_series_2017_2018_S2A_SP.tif'
Datefile = 'NDVI_BXII_time_series_images_list.csv'
# open GeoTIFF dataset as an object
src = rasterio.open(path + NDVIfile)
# print meta data (selection)
for key in src.meta:
print(key,': ',src.meta[key])
print('\nbounds',src.bounds)
print('wkt',src.crs.wkt)
# compute imshow extent
extent = [ src.bounds.left , src.bounds.right , src.bounds.bottom , src.bounds.top ]
print('extent',extent)
# Read the image data
NDVI = src.read()
print('\nNDVI shape',NDVI.shape)
driver : GTiff dtype : float64 nodata : None width : 866 height : 800 count : 33 crs : EPSG:32630 transform : | 20.00, 0.00, 214920.00| | 0.00,-20.00, 4106200.00| | 0.00, 0.00, 1.00| bounds BoundingBox(left=214920.0, bottom=4090200.0, right=232240.0, top=4106200.0) wkt PROJCS["WGS 84 / UTM zone 30N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-3],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32630"]] extent [214920.0, 232240.0, 4090200.0, 4106200.0] NDVI shape (33, 800, 866)
# Read dates using Pandas
datetable = pd.read_csv(path + Datefile)
# Preview the first 5 lines of the loaded dates
print("The date table .head() look like this: \n",datetable.head())
# Preview the data types
print(datetable.dtypes)
# Convert to 1-D numpy arrays
ts = np.array(datetable.Day)
tsdatestrings = np.array(datetable.ISOdate)
tsdoy = np.array(datetable.Doy)
print('ts',ts)
print('tsdatestrings',tsdatestrings)
print('tsdoy',tsdoy)
datetable
The date table .head() look like this:
Band ISOdate Day Doy
0 1 2017-04-12 1 102
1 2 2017-04-22 11 112
2 3 2017-05-02 21 122
3 4 2017-05-12 31 132
4 5 2017-05-22 41 142
Band int64
ISOdate object
Day int64
Doy int64
dtype: object
ts [ 1 11 21 31 41 51 61 71 81 91 101 111 131 141 151 161 171 181
191 201 211 221 231 241 261 271 291 301 311 321 341 351 371]
tsdatestrings [' 2017-04-12' ' 2017-04-22' ' 2017-05-02' ' 2017-05-12' ' 2017-05-22'
' 2017-06-01' ' 2017-06-11' ' 2017-06-21' ' 2017-07-01' ' 2017-07-11'
' 2017-07-21' ' 2017-07-31' ' 2017-08-20' ' 2017-08-30' ' 2017-09-09'
' 2017-09-19' ' 2017-09-29' ' 2017-10-09' ' 2017-10-19' ' 2017-10-29'
' 2017-11-08' ' 2017-11-18' ' 2017-11-28' ' 2017-12-08' ' 2017-12-28'
' 2018-01-07' ' 2018-01-27' ' 2018-02-06' ' 2018-02-16' ' 2018-02-26'
' 2018-03-18' ' 2018-03-28' ' 2018-04-17']
tsdoy [102 112 122 132 142 152 162 172 182 192 202 212 232 242 252 262 272 282
292 302 312 322 332 342 362 7 27 37 47 57 77 87 107]
| Band | ISOdate | Day | Doy | |
|---|---|---|---|---|
| 0 | 1 | 2017-04-12 | 1 | 102 |
| 1 | 2 | 2017-04-22 | 11 | 112 |
| 2 | 3 | 2017-05-02 | 21 | 122 |
| 3 | 4 | 2017-05-12 | 31 | 132 |
| 4 | 5 | 2017-05-22 | 41 | 142 |
| 5 | 6 | 2017-06-01 | 51 | 152 |
| 6 | 7 | 2017-06-11 | 61 | 162 |
| 7 | 8 | 2017-06-21 | 71 | 172 |
| 8 | 9 | 2017-07-01 | 81 | 182 |
| 9 | 10 | 2017-07-11 | 91 | 192 |
| 10 | 11 | 2017-07-21 | 101 | 202 |
| 11 | 12 | 2017-07-31 | 111 | 212 |
| 12 | 13 | 2017-08-20 | 131 | 232 |
| 13 | 14 | 2017-08-30 | 141 | 242 |
| 14 | 15 | 2017-09-09 | 151 | 252 |
| 15 | 16 | 2017-09-19 | 161 | 262 |
| 16 | 17 | 2017-09-29 | 171 | 272 |
| 17 | 18 | 2017-10-09 | 181 | 282 |
| 18 | 19 | 2017-10-19 | 191 | 292 |
| 19 | 20 | 2017-10-29 | 201 | 302 |
| 20 | 21 | 2017-11-08 | 211 | 312 |
| 21 | 22 | 2017-11-18 | 221 | 322 |
| 22 | 23 | 2017-11-28 | 231 | 332 |
| 23 | 24 | 2017-12-08 | 241 | 342 |
| 24 | 25 | 2017-12-28 | 261 | 362 |
| 25 | 26 | 2018-01-07 | 271 | 7 |
| 26 | 27 | 2018-01-27 | 291 | 27 |
| 27 | 28 | 2018-02-06 | 301 | 37 |
| 28 | 29 | 2018-02-16 | 311 | 47 |
| 29 | 30 | 2018-02-26 | 321 | 57 |
| 30 | 31 | 2018-03-18 | 341 | 77 |
| 31 | 32 | 2018-03-28 | 351 | 87 |
| 32 | 33 | 2018-04-17 | 371 | 107 |
To plot the NDVI for the first day we are going to use matplotlib's imshow method to make a pseudo-color map. A good colormap for the NDVI is essential! It is up to you to come up with a good colormap for the NDVI, taking into account the lessons learned from Exercise 1 and 2. You may also use your answer from Exercise 4 here (this is a repetition).
Start by formulating the requirements for the color map. What kind of color map should we select? What is the range? Do we use a linear map, is it continuous or discrete, do we use two tones? What kind of normalization are we going to use? In any case, the default color map is certainly not a good choice. Why?
Have a look for instance at the following article https://publiclab.org/notes/cfastie/08-26-2014/new-ndvi-colormap to get some inspiration.
In the section below set the cmap variable with the chosen colormap.
In Python you can read the text file or mat file and create the colormap using matplotlib.colors.ListedColormap. See also Exercise 1 and 2.
cmapTable = scipy.io.loadmat('file.mat')['colors']
cmap = colors.ListedColormap(data/255)
You can then declare this new colormap in the imshow function by declaring the input as cmap = cmap
# SELECT APPROPRIATE COLORMAP FOR NDVI
cmap = 'RdYlGn' # default colormap for matplotlib, change this for something better!!
# END COLORMAP SELECTION
clim = [-1 , 1] # NDVI range
To plot the NDVI for the first day we are going to use matplotlib's imshow method to make a pseudo-color map.
Use the colormap that you selected above using the cmap = cmap optional parameter in imshow, together with 'clim' and 'extent'.
The optional extent parameter is used to set the proper x- and y-axis coordinates. It is computed from the image bounds, while converting meter units to km.
# Put your code here ((you may borrow it from Exercise 4)
extentKiloMeter = np.array(extent)/1000
# borrowed from exercise 4
%matplotlib inline
band = 0
plt.figure(figsize=(8, 8))
fig = plt.imshow(NDVI[band, :, :], extent=extentKiloMeter, cmap=cmap, clim=clim)
plt.title('NDVI 2017-04-12')
plt.title('NDVI ' + tsdatestrings[band] + ' DAY ' + str(int(ts[band])))
plt.ylabel('North [km]')
plt.xlabel('East [km]')
plt.colorbar(label='NDVI [-]', fraction=0.0417, pad=0.06)
plt.show()
A video is created showing the (unfiltered) NDVI data. We use imshow to display matrix NDVI(i,:,:) as an image for day i, using the colormap of your choice (see previous section).
The figure statement is called only once, imshow is called from within the loop over days, and the output is appended to the list of frames.
There are two ways to display the animation from a notebook,
You can also save the animation as either a mpeg file or an animated gif (both are commented out at the moment). The code is borrowed from Exercise 4.
# Put your code here, you may borrow it from Exercise 4
# // as is copy paste from exercise 4
# %matplotlib qt # for using Qt
# %matplotlib notebook # for Qt inline in notebook
# %matplotlib inline # for inline html5
fig, ax=plt.subplots(figsize=(10,8))
# frames is a list of lists, each row is a list of artists to draw in the current frame;
# here we are just animating one artist, the image, in each frame
frames = []
ndates = NDVI.shape[0]
for i in range(ndates):
titlestr = 'NDVI ' + tsdatestrings[i] + ' DAY ' + str(int(ts[i])) # title = Add your title here
im = ax.imshow(NDVI[i, :, :], extent=extentKiloMeter, aspect='equal', cmap=cmap, clim=clim, animated=True) # im = Add your images here
title = ax.text(0.5, 1.05, titlestr, size=plt.rcParams["axes.titlesize"], ha="center", transform=ax.transAxes, animated=True)
frames.append([im, title])
plt.ylabel('North [km]')
plt.xlabel('East [km]')
fig.colorbar(im, ax=ax, label='NDVI [-]', fraction=0.0417, pad=0.06)
ani = animation.ArtistAnimation(fig, frames, interval=300, blit=False, repeat_delay=500)
FFwriter = animation.FFMpegWriter(fps=3)
ani.save('NDVI_raw.mp4', writer = FFwriter)
The colormap used in the above animation and the plot is a diverging colormap from red to green. It has been chosen to portray that the greener the pixel is, the more NDVI value it has and vice versa with respect to the reds in the map.
Inputs for HANTS are the times ts (days) and the NDVI time series for one pixel. For the other input parameters see the help that comes with HANTS.
help(hants)
Help on function hants in module hants:
hants(ni, nb, nf, y, ts, HiLo, low, high, fet, dod, delta)
Harmonic ANalysis of Time Series (HANTS)
Parameters
----------
ni : int
Number of images (total number of actual samples of the time series)
nb : int
Length of the base period, measured in virtual samples (days, dekads, months, etc.)
nf: int
Number of frequencies to be considered above the zero frequency
y: ndarray
1D array of size size `ni` with input sample values (e.g. NDVI values)
ts: ndarray
1D array of size `ni` of time sample indicators (indicates virtual
sample number relative to the base period); numbers in array `ts`
maybe greater than `nb`
HiLo: str {'Hi', 'Lo'}
2-character string indicating rejection of high 'Hi' or low 'Lo' outliers
low: float
valid range minimum
high: float
valid range maximum (values outside the valid range are rejeced right away)
fet: float
fit error tolerance (points deviating more than fet from curve fit are rejected)
dod : int
degree of overdeterminedness (iteration stops if number of points reaches
the minimum required for curve fitting, plus dod). This is a safety measure
delta : float
small positive number (e.g. 0.1) to suppress high amplitudes
Returns
-------
amp : ndarray
1D array of length `nf`+` with Amplitude for frequencies 0...nf
phi : ndarray
1D array of length `nf'+1 with Phase in degrees for frequencies 0...nf
yr : ndarray
1D array of length `ni` with smoothed values (e.g. NDVI)
outliers : ndarray
1D Boolean array of length `ni` with flagged outliers
Notes
-----
This function applies the Harmonic ANalysis of Time Series (HANTS)
algorithm originally developed by the Netherlands Aerospace Centre (NLR)
(http://www.nlr.org/space/earth-observation/).
This python implementation was based on two previous implementations
available at the following links:
https://codereview.stackexchange.com/questions/71489/harmonic-analysis-of-time-series-applied-to-arrays
http://nl.mathworks.com/matlabcentral/fileexchange/38841-matlab-implementation-of-harmonic-analysis-of-time-series--hants-
Original Author for Python Version: Espinoza-Dávalos, G. E., Bastiaanssen, W. G. M., Bett, B., & Cai, X. (2017).
A Python Implementation of the Harmonic ANalysis of Time Series (HANTS) Algorithm for Geospatial Data.
http://doi.org/10.5281/zenodo.820623
Edited by: Ullas Rajvanshi and Hans van der Marel
Based on the help information define your HANTS input parameters below. Some are preset, others are left open.
nb= 371 # length of the base period (usually one year) measured in timeunits of ts (number of days)
nf= 1 # number of frequencies to be considered above the zero frequency ( 1 ... 6)
fet= 0.05 # fit error tolerance (points deviating more than fet from curve fit are rejected, 0.05 ... 0.3)
HiLo='Lo' # 2-character string indicating rejection of high or low outliers ('Lo' or 'Hi')
low=-1 # valid range minimum (values below are rejected right away)
high=+1 # valid range maximum (values above are rejected right away)
dod=1 # degree of overdeterminedness (iteration stops if number of points reaches the minimum required)
delta=0.1 # small positive number (e.g. 0.1) to suppress high amplitudes
Set the line and sample number of the pixel for which to apply HANTS.
We will show you later how to find the line and sample numbers for specific pixel areas, e.g. one of the field in the training set.
row = 348
col = 299
Do HANTS, first we have to extract the pixel data into the array datain. We use the Matlab function squeeze to do this, because NDVI is stored as a three dimensional array and we have to collect all values that belong to a single pixel. We use squeeze because NDVI[:,row,col] does not return an array, but a list, and squeeze will convert this list to an array.
HANTS Outputs are:
amp = array of amplitudes, first element is the average of the curve
phi = array of phases, first element is zero
dataout = array holding reconstructed time series
outliers = boolean array with the outliers
ndates = NDVI.shape[0]
datain = np.squeeze(NDVI[:, row, col])
amp, phi, dataout, outliers = hants(ndates, nb, nf, datain, ts, HiLo, low, high, fet, dod, delta)
print('amp',amp)
print('phi',phi)
print('dataout',dataout)
print('outliers',outliers)
numoutliers=np.sum(outliers)
print('numoutliers',numoutliers)
amp [0.49355339 0.40122089] phi [ 0. 51.22445148] dataout [0.75008759 0.79841392 0.83801707 0.86776384 0.88680308 0.89459001 0.8909018 0.87584399 0.84984744 0.81365601 0.76830527 0.71509287 0.59135488 0.52436988 0.45650312 0.3896965 0.3258616 0.26682499 0.2142759 0.16971797 0.13442617 0.10941031 0.09538619 0.09275511 0.121645 0.15233934 0.24186053 0.29812585 0.35998307 0.42566224 0.5609135 0.62661556 0.74482646] outliers [ True True True True True True False False False False True True True True True True True True True True False False True False True True False True True True True True True] numoutliers 25
%matplotlib inline
# making a function to replicate the plotting
def readymade_HANTS_plot_maker(ts, datain, outliers, dataout, row, col, nb, nf, fet):
numoutliers=np.sum(outliers)
plt.figure(figsize=(10,6))
plt.plot(ts, datain, 'b*', label='Original NDVI')
plt.plot(ts[outliers], datain[outliers], 'rx', markersize=10, markeredgewidth=3, label='Outlier NDVI')
plt.plot(ts, dataout, 'r', label='Smoothed with HANTS')
plt.xlim([1, 365])
#plt.ylim([0, 1])
plt.xlabel('Days from 12 April 2017')
plt.ylabel('NDVI [-]')
plt.title('NDVI - Sentinel 2 (row {}, column {})'.format(row, col));
plt.text(6,0.77,'nb={}\nnf={}\nfet={}\noutliers={}'.format(nb, nf, fet, numoutliers))
plt.legend()
plt.show()
readymade_HANTS_plot_maker(ts, datain, outliers, dataout, row, col, nb, nf, fet)
this will help us to iterate through different values of frequencies and days
# check hants with different frequencies
def HANTS_and_plottting(ndates, nb, nf, datain, ts, HiLo, low, high, fet, dod, delta, row,col):
ndates = NDVI.shape[0]
datain = np.squeeze(NDVI[:, row, col])
amp, phi, dataout, outliers = hants(ndates, nb, nf, datain, ts, HiLo, low, high, fet, dod, delta)
numoutliers=np.sum(outliers)
readymade_HANTS_plot_maker(ts, datain, outliers, dataout, row, col, nb, nf, fet)
HANTS_and_plottting(ndates, nb, 3, datain, ts, HiLo, low, high, 10/100, dod, delta, row, col)
HANTS_and_plottting(ndates, nb, 3, datain, ts, HiLo, low, high, 10/100, dod, delta, 200,150)
# for i in range(1, 7, 1):
# # changing frequency for everything remaining same
# for j in range(5, 31, 5) :
# HANTS_and_plottting(ndates, nb, i, datain, ts, HiLo, low, high, j/100, dod, delta)
# print(f"Changing frequency to {i+1}\n")
# # HANTS_and_plottting(ndates, nb, nf, datain, ts, HiLo, low, high, , dod, delta)a
Use the above code as a template to select a suitable set of input parameters. You can try different values for the number of frequencies (nf) , base period (nb) and fit error tolerance.
You should at least
In the report, explain your choice for the base period (you don't have to show plots, there is only one good value for the base period), and include the plots for the chosen parameters.
# Add your plots here
nb=371
nf=3
fet=0.1
Base period: It is chosen as 371 since the 33 NDVI images are over the period of 371 days.
nf: 3
fet: 0.1
From the above plots, it is decided to take nf value as 3 and the tolerance set at 0.1, since at nf=4 we start seeing that the curve almost fits the highest values whereas on increasing the nf value further it is observed that the curve can be called as over-fitting the values.
Apply the HANTS algorithm for all pixels, and save the amplitude, phase and reconstucted NDVI.
print('Smoothing and filling the missing data, please be patient ... ')
# Initialize new arrays to collect the HANTS output
ni, ny, nx = NDVI.shape
NDVI_HANTS = np.zeros((ni, ny, nx))
amp = np.zeros((nf + 1, ny, nx))
phi = np.zeros((nf + 1, ny, nx))
# Call HANTS in a double loop, line by line (SEE THE WARNING)
#for row in range(0, 10): # use this for smaller sample
for row in range(0, ny):
if ( row % 25 ) == 0: print("Doing line {} out of {}".format(row, ny))
for col in range(0, nx):
data = NDVI[:, row, col]
if np.isnan(data).sum() != ni:
data[np.isnan(data)] = low - 1.0
amp_, phi_, yr_, outlier_ = hants(ndates, nb, nf, data, ts, HiLo, low, high, fet, dod, delta)
amp[:, row, col], phi[:, row, col], NDVI_HANTS[:, row, col] = amp_, phi_, yr_
print("Done with smoothing and filling")
Smoothing and filling the missing data, please be patient ... Doing line 0 out of 800 Doing line 25 out of 800 Doing line 50 out of 800 Doing line 75 out of 800 Doing line 100 out of 800 Doing line 125 out of 800 Doing line 150 out of 800 Doing line 175 out of 800 Doing line 200 out of 800 Doing line 225 out of 800 Doing line 250 out of 800 Doing line 275 out of 800 Doing line 300 out of 800 Doing line 325 out of 800 Doing line 350 out of 800 Doing line 375 out of 800 Doing line 400 out of 800 Doing line 425 out of 800 Doing line 450 out of 800 Doing line 475 out of 800 Doing line 500 out of 800 Doing line 525 out of 800 Doing line 550 out of 800 Doing line 575 out of 800 Doing line 600 out of 800 Doing line 625 out of 800 Doing line 650 out of 800 Doing line 675 out of 800 Doing line 700 out of 800 Doing line 725 out of 800 Doing line 750 out of 800 Doing line 775 out of 800 Done with smoothing and filling
# # Save important variables
import pickle
# f = open('NDVIr.pckl', 'wb')
# pickle.dump([ts, tsdatestrings, extent, NDVI_HANTS, amp, phi], f)
# f.close()
If the variables have to be reloaded you can use the following code.
# ReLoad the variables
f = open('NDVIr.pckl', 'rb')
ts, tsdatestrings, extent, NDVI_HANTS, amp, phi = pickle.load(f)
f.close()
Three new GeoTIFF files are created, with
These files can easily be imported by other software such as QGIS.
# Save reconstructed NDVI as GeoTiff file
profile = src.profile
with rasterio.open('NDVI_reconstructed.tif', 'w', **profile) as dst:
for id, layer in enumerate(NDVI_HANTS, 1):
dst.write_band(id, NDVI_HANTS[id - 1, :, :])
# Save Amplitude as GeoTiff file
profile = src.profile
profile['count']=len(amp)
with rasterio.open('NDVI_amplitude.tif', 'w', **profile) as dst:
for id, layer in enumerate(amp, 1):
dst.write_band(id, amp[id - 1, :, :])
# Save Phase as GeoTiff file
profile = src.profile
profile['count']=len(phi)
with rasterio.open('NDVI_phase.tif', 'w', **profile) as dst:
for id, layer in enumerate(phi, 1):
dst.write_band(id, phi[id - 1, :, :])
Create a video showing the reconstructed NDVI data.
# Your code here
# create a function to return a list of image and title
def dayMapper(no_day):
reconst_date_ = tsdatestrings[no_day]
img_local_ax = ax.imshow( NDVI_HANTS[no_day, :, :], cmap=cmap, clim=clim, animated=True, aspect="equal" )
img_title = ax.text(150, -20, f"NDVI Reconstructed with HANTS on {reconst_date_}")
return [img_local_ax, img_title]
fig, ax = plt.subplots()
ndates = NDVI_HANTS.shape[0]
# make frames
frames = [dayMapper(i) for i in range(1, ndates)]
# generate animation object
anim = animation.ArtistAnimation( fig, frames, interval=300, blit=False, repeat_delay=500 )
# save resultsF
fig.colorbar(frames[0][0], ax=ax, fraction=0.0417, pad=0.06)
fig.set_facecolor("white")
anim.save('NDVI_reconstructed_animation.mp4', writer = FFwriter)
Instead of an animation we can also create a lot of subplots using the code below.
%matplotlib inline
ndates = NDVI_HANTS.shape[0]
ncols = 5
nrows = int(np.ceil(ndates/ncols))
fig, axs = plt.subplots(nrows, ncols, figsize=(21,27), tight_layout=True)
axs= axs.flat
[axi.set_axis_off() for axi in axs]
for i in range(ndates):
titlestr = 'Reconstructed NDVI ' + tsdatestrings[i] + ' DAY ' + str(int(ts[i]))
axs[i].imshow(NDVI_HANTS[i, :, :], extent=extent, aspect='equal', cmap=cmap, clim=clim)
axs[i].set_xticks([])
axs[i].set_yticks([])
axs[i].set_title(titlestr)
plt.show()
Make a figure with in the subplots a pseudo color image of the amplitude and phase for each frequency. Please think carefully about the colormap to use.
Which colormap do you use for the amplitude. Should the colormap for the phase be the same, or different. Please explain why.
# Your code here
plt.figure(figsize=(15,15))
# first plot
plt.subplot(1, 2, 1)
amplitude_cmap = 'PuBuGn'#'YlGnBu'
plt.imshow(amp[1, :, :], cmap=amplitude_cmap, label="Amplitude", extent=extentKiloMeter)
plt.colorbar(fraction=0.042)
plt.ylabel("North [km]")
plt.xlabel("East [km]")
plt.title("Amplitude of reconstructed time series")
# 2nd plot
plt.subplot(1, 2, 2)
phase_cmap = 'YlGnBu'#'PuBuGn'
plt.imshow(phi[1, :, :], cmap=phase_cmap, extent=extentKiloMeter)
plt.colorbar(fraction=0.042, label="Phase angle [Degrees]")
plt.title("Phase angle of reconstructed time series")
plt.ylabel("North [km]")
plt.xlabel("East [km]")
plt.show()
28 November 2017 was significantly cloudy, so we will use that image, from the dataframe we can determine that the index for this image is 22
plt.figure(figsize=(15,15))
# first plot
plt.subplot(1, 2, 1)
plt.imshow(NDVI[22, :, :], cmap=cmap, extent=extentKiloMeter)
plt.colorbar(fraction=0.042)
plt.ylabel("North [km]")
plt.xlabel("East [km]")
plt.title("Unfiltered NDVI 2017-11-28")
# 2nd plot
plt.subplot(1, 2, 2)
plt.imshow(NDVI_HANTS[22, :, :], cmap=cmap, extent=extentKiloMeter)
plt.colorbar(fraction=0.042, label="NDVI")
plt.ylabel("North [km]")
plt.xlabel("East [km]")
plt.title("Reconstructed NDVI 2017-11-28")
plt.show()
The training dataset is a shapefile with information on the actual crop that is growing in the fields. The original training dataset was in UTM zone 29N. We used QGIS to convert this to UTM zone 30N, which was is also the zone that is used for the NDVI data.
First we must read the shape file with training data
training = shapefile.Reader("Training_set/Training_2017_UTM30N_WGS84_SP.shp")
# Your code here
training = shapefile.Reader("./input_files/Training_set/Training_2017_UTM30N_WGS84_SP.shp")
training
<shapefile.Reader at 0x7f46f4cefa60>
The training dataset consists of 15 polygons. We will create an plot with the amplitude data and overlay the polygon boundaries.
plt.figure(figsize=(12,8))
for shape in training.shapeRecords():
x = [i[0]/1000 for i in shape.shape.points[:]]
y = [i[1]/1000 for i in shape.shape.points[:]]
plt.plot(x,y, color='yellow', linewidth=2)
greys = mpl.cm.get_cmap('Greys_r')
plt.imshow(amp[0, :, :], extent=extentKiloMeter, cmap=greys)
plt.axis('image')
plt.colorbar(label='Amplitude [-]')
plt.ylabel('North [km]')
plt.xlabel('East [km]')
plt.title('NDVI Training set')
plt.show()
The amplitude obtained from the reconstruction of the images was loaded into QGIS and the amplitudes stored in bands 1,3,4 were mapped to the R,G,B bands on the visible image.
The image finally looks like the following:

[End of the assignment]